This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present invention. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present invention. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.
Viscoacoustic and Viscoelastic Wavefield Modeling
In the time-domain, intrinsic attenuation (absorption and dispersion) is formulated with convolution operators between strain and relaxation functions. Such convolutions are computationally impractical for large-scale wave propagation simulations using time-marching methods. Incorporation of realistic attenuation into time-domain computations was first achieved using Padé approximants by Day and Minster (Day and Minster, 1984). Later, Padé approximants evolved into sophisticated viscoacoustic and viscoelastic rheological models in rational forms, such as the generalized-Maxwell (GMB) and Standard-Linear-Solid (SLS) models. Attenuation modeling methods in the time domain are based on either the GMB or SLS formulations given by Emmerich and Korn (1987) and Carcione et al. (1988), respectively, and use rheological models based on relaxation mechanisms.
A relaxation mechanism is the unit of the time-domain attenuation model representing viscous effects in a narrow frequency band. Multiple relaxation mechanisms are combined to model attenuation over a desired frequency band, introducing additional state variables and partial differential equations (PDE) to the acoustic and elastic forward wave equations. A considerable amount of computational time and memory are consumed by these additional variables and equations during forward wave simulations and to an even greater degree during adjoint simulations.
Only simplified versions of approximate SLS formulations have been incorporated into FWI to infer a frequency-invariant attenuation at a reduced computational cost, such as Robertsson et al. 1994, Charara et al. 2000, Hestholm et al. 2006, and Royle 2011. Not only does the accuracy of such approximations deteriorate with decreasing quality factor values (stronger attenuation), but also they cannot be used to model quality factors that vary with frequency. Another approach is to limit the number of relaxation mechanisms to one or two in order to minimize the complexity of the SLS model, which results in poor representation of the targeted quality factor over the seismic frequency band of interest. For example, Charara (Charara 2006) suggests using only two relaxation mechanisms for the SLS model to reduce the computational cost in computing gradients. Their formulation requires computing spatial derivatives of memory variables in the adjoint simulations, which limits them in the total number of relaxation mechanisms that can be used in their approximate SLS model. Also, Bai (Bai et al., 2012) presented a method for compensating attenuation effects for FWI using a single mechanism SLS model with the assumption of a frequency-invariant quality factor.
While only SLS attenuation models have been incorporated into time-domain viscoacoustic and viscoelastic FWI techniques, both SLS and GMB models have been used for forward viscoacoustic and viscoelastic modeling of waves. For example, Käser (Käser et al., 2007) applied GMB for modeling viscoelastic waves with a discontinuous Galerkin method.
A common measure of attenuation is the quality factor, a dimensionless quantity that defines the frequency dependence of the acoustic or elastic moduli. The quality factor itself can be frequency dependent, especially for fluid-bearing rocks such as hydrocarbons, and is typically assumed to be frequency invariant for dry rocks (Muller et al., 2010; Quintal, 2012). Therefore, it is useful to infer the frequency dependence of the quality factor to identify fluid bearing rocks such as reservoirs.
The SLS attenuation model represents the numerical inverse of the quality factor Q, also called the loss factor, by the following relation (Carcione et al., 1988):
                                          Q                          -              1                                ⁡                      (                          x              ,              ω                        )                          =                                            ∑                              l                =                1                            L                        ⁢                                                  ⁢                                          ω                ⁡                                  (                                                                                    τ                                                  ɛ                          ⁢                                                                                                          ⁢                          l                                                                    ⁡                                              (                        x                        )                                                              -                                                                  τ                                                  σ                          ⁢                                                                                                          ⁢                          l                                                                    ⁡                                              (                        x                        )                                                                              )                                                            1                +                                                                            ω                      2                                        ⁡                                          (                                                                        τ                                                      σ                            ⁢                                                                                                                  ⁢                            l                                                                          ⁡                                                  (                          x                          )                                                                    )                                                        2                                                                                        ∑                              l                =                1                            L                        ⁢                                                  ⁢                                          1                +                                                      ω                    2                                    ⁢                                                            τ                                              ɛ                        ⁢                                                                                                  ⁢                        l                                                              ⁡                                          (                      x                      )                                                        ⁢                                                            τ                                              σ                        ⁢                                                                                                  ⁢                        l                                                              ⁡                                          (                      x                      )                                                                                                  1                +                                                                            ω                      2                                        ⁡                                          (                                                                        τ                                                      σ                            ⁢                                                                                                                  ⁢                            l                                                                          ⁡                                                  (                          x                          )                                                                    )                                                        2                                                                                        (        1        )                            where        Q=quality factor,        τϵl=strain relaxation time of mechanism l in SLS model,        τσl=stress relaxation time of mechanism l in SLS model,        x=spatial coordinate,        ω=frequency,        L=number of relaxation mechanisms used in the SLS model.        Conceptually, the quality factor Q represents the ratio of stored to dissipated energy in a medium. The strain and stress relaxation times are determined to best fit the desired quality factor distribution over the frequency band.        
A conventional first-order form of the linear viscoacoustic wave equations for simulating waves in attenuating acoustic media (JafarGandomi et al. 2007) is:
                                                                        ∂                p                                            ∂                t                                      +                          κ              ⁢                              ∇                                  ·                  v                                                      +                                          ∑                                  l                  =                  1                                L                            ⁢                                                          ⁢                                                ϕ                  l                                ⁢                                  m                  l                                                              =                      s            p                          ,                                  ⁢                                                            ∂                v                                            ∂                t                                      +                                          1                ρ                            ⁢                              ∇                p                                              =                      s            v                          ,                                  ⁢                                                            ∂                                  m                  l                                                            ∂                t                                      +                                          κα                l                            ⁢                              ∇                                  ·                  v                                                      +                                          ϕ                l                            ⁢                              m                l                                              =          0                ,                            (        2        )                            with appropriate initial and boundary conditions for pressure p, velocity v, and memory variables ml. Note that        ∇=divergence operator,        κ=unrelaxed bulk modulus (limω→∞κ(ω)→κ),        ρ=mass density,        v=velocity (v={vx vy vz}T in 3D space),        p=pressure,        ml=memory variable for mechanism l,        sp=pressure source,        sv=velocity source,        
      ϕ    l    =                    1                  τ                      σ            ⁢                                                  ⁢            l                              ⁢                          ⁢      and      ⁢                          ⁢              α        l              =          (              1        -                              τ                          σ              ⁢                                                          ⁢              l                                            τ                          ɛ              ⁢                                                          ⁢              l                                          )                      where relaxation parameters τϵl and τσl may be determined by Equation (1) for a given quality factor profile.        Note that continuous scalar variables are denoted by italicized characters and vector and matrices are denoted by bold non-italicized characters throughout this document.        
FWI methods based on computing gradients of an objective function with respect to the parameters are often efficiently implemented by using adjoint methods, which have been proved to outperform other relevant methods, such as direct sensitivity analyses, finite differences or complex variable methods. One can find a detailed discussion about adjoint methods in Thevenin et al., 2008.
The continuous adjoint of the conventional viscoacoustic system (Equations (2)) is
                                                                        ∂                                  p                  _                                                            ∂                t                                      +                          ∇                              ·                                  (                                                            1                      ρ                                        ⁢                                          v                      _                                                        )                                                              =                                    ∂              ℱ                                      ∂              p                                      ,                                  ⁢                                                            ∂                                  v                  _                                                            ∂                t                                      +                          ∇                              (                                  κ                  ⁢                                                                          ⁢                                      p                    _                                                  )                                      +                                          ∑                                  l                  =                  1                                L                            ⁢                                                          ⁢                              ∇                                  (                                                            κα                      l                                        ⁢                                                                  m                        _                                            l                                                        )                                                              =                                    ∂              ℱ                                      ∂              v                                      ,                                  ⁢                                                            ∂                                                      m                    _                                    l                                                            ∂                t                                      +                                          ϕ                l                            ⁢                              p                _                                      +                                          ϕ                l                            ⁢                                                m                  _                                l                                              =          0                ,                            (        3        )                            where        p=adjoint pressure,        v=adjoint velocity,        ml=adjoint memory variable for mechanism l, and∂/∂p and ∂/∂v are derivatives of the objective function  with respect to the pressure and velocity respectively. Although the objective function and its derivative with respect to the pressure and velocity will be explicitly defined in the Viscoacoustic and Viscoelastic Full Wavefield Inversion Section, ∂/∂p and ∂/∂v can simply be considered as sources to the adjoint equations.        
Note that the derivation of the adjoint equations involves using integration by parts, which introduces both spatial and temporal boundary terms into the adjoint equations (which are evaluated at the spatial and temporal boundaries respectively). These adjoint spatial boundary terms are not included in Equations (3), but need to be included (i.e. in code) in order to correctly compute the unique solutions for the adjoint variables, and the temporal boundary terms are zero and thus drop out.
The spatial derivatives of memory variables ∇καlml in the adjoint equations (Equations (3)) lead to a significant amount of computational cost. FIG. 5 shows the cost of adjoint computations relative to forward computations for second through twelfth order (spatially) accurate finite-difference (FD) time-domain methods. The relative cost is based on comparing the number of floating-point operations (flops). For three relaxation mechanisms and an eighth-order finite-difference method, the cost of solving the adjoint equations is 2.25 times more expensive than the cost of solving the forward equations (Equations (2)). Also note that as the spatial order of the finite-difference method and the number of relaxation mechanisms increase, the relative cost of solving the adjoint equations increases. Furthermore, it is expected that for large-scale distributed-memory parallel processing computations, the cost of solving the adjoint will increase due to the additional communication of the memory variables.
The linear viscoacoustic equations can be extended to the following linear viscoelastic equations:
                                                                        ∂                σ                                            ∂                t                                      +                          C              ⁢                                              T                            ⁢              v                        +                                          ∑                                  l                  =                  1                                L                            ⁢                                                          ⁢                                                ϕ                  l                                ⁢                                  m                  l                                                              =                      s            σ                          ,                                  ⁢                                                            ∂                v                                            ∂                t                                      +                                          1                ρ                            ⁢                            ⁢              σ                                =                      s            v                          ,                                  ⁢                                                            ∂                                  m                  l                                                            ∂                t                                      +                                          D                l                            ⁢                                              T                            ⁢              v                        +                                          ϕ                l                            ⁢                              m                l                                              =          0                ,                            (        4        )                            where in the 3D Cartesian coordinate system,        σ=stress, σ={σxx σyy σzz σxz σyz σxy}T,        ml=memory variable for mechanism l,        ml={ml,xx ml,yy ml,zz ml,xz ml,yz ml,xy}T,        
      C    =          [                                                  λ              +                              2                ⁢                μ                                                          λ                                λ                                0                                0                                0                                                λ                                              λ              +                              2                ⁢                μ                                                          λ                                0                                0                                0                                                λ                                λ                                              λ              +                              2                ⁢                μ                                                          0                                0                                0                                                0                                0                                0                                μ                                0                                0                                                0                                0                                0                                0                                μ                                0                                                0                                0                                0                                0                                0                                μ                              ]        ,                which is the elastic constitutive relationship for the isotropic unrelaxed system in terms of Lamé constants λ and μ,        
            D      l        =          [                                                                  λα                l                λ                            +                              2                ⁢                                  μα                  l                  μ                                                                                        λα              l              λ                                                          λα              l              λ                                            0                                0                                0                                                              λα              l              λ                                                                          λα                l                λ                            +                              2                ⁢                                  μα                  l                  μ                                                                                        λα              l              λ                                            0                                0                                0                                                              λα              l              λ                                                          λα              l              λ                                                                          λα                l                λ                            +                              2                ⁢                                  μα                  l                  μ                                                                          0                                0                                0                                                0                                0                                0                                              μα              l              μ                                            0                                0                                                0                                0                                0                                0                                              μα              l              μ                                            0                                                0                                0                                0                                0                                0                                              μα              l              μ                                          ]        ,                which is the constitutive relationship for the memory system (for the SLS formulation        
            α      l              λ        ,        μ              =          (              1        -                              τ                          σ              ⁢                                                          ⁢              l                                            τ                          ɛ              ⁢                                                          ⁢              l                                      λ              ,              μ                                          )        ,                and τϵlλ and τϵlμ are computed from the compressional and shear wave quality factors QP and QS using Equation (1)),        T=strain operator,        
          T    =                    [                                                            ∂                x                                                    0                                      0                                                      ∂                z                                                    0                                                      ∂                y                                                                        0                                                      ∂                y                                                    0                                      0                                                      ∂                z                                                                    ∂                x                                                                        0                                      0                                                      ∂                z                                                                    ∂                x                                                                    ∂                y                                                    0                                      ]            T        .                  The continuous adjoint of the conventional viscoelastic system (Equations (4)) is        
                                                                        ∂                                  σ                  _                                                            ∂                t                                      +                                                          T                            ⁢                              (                                                      1                    ρ                                    ⁢                                      v                    _                                                  )                                              =                                    ∂              ℱ                                      ∂              σ                                      ,                                  ⁢                                                            ∂                                  v                  _                                                            ∂                t                                      +                                        ⁢                              (                                  C                  ⁢                                                                          ⁢                                      σ                    _                                                  )                                      +                                          ∑                                  l                  =                  1                                L                            ⁢                                                          ⁢                                              ⁢                                  (                                                            D                      l                                        ⁢                                                                  m                        _                                            l                                                        )                                                              =                                    ∂              ℱ                                      ∂              v                                      ,                                  ⁢                                                            ∂                                                      m                    _                                    l                                                            ∂                t                                      +                                          ϕ                l                            ⁢                              σ                _                                      +                                          ϕ                l                            ⁢                                                m                  _                                l                                              =          0                ,                            (        5        )                            where in the 3D Cartesian coordinate system,        σ=adjoint stress, σ={σxx σyy σzz σxz σyz σxy}T,        ml=adjoint memory field for mechanism l,        ml={ml,xx ml,yy ml,zz ml,xz ml,yz ml,xy}T, and ∂/∂σ and ∂/∂v are derivatives of the objective function  with respect to the stress and velocity respectively.        
Note that the spatial boundary terms are not included in Equations (4). The temporal boundary terms are zero and thus drop out.
As for the viscoacoustic case, the adjoint equations for the viscoelastic case also contain spatial derivatives of memory variables, which will add a large computational cost to viscoelastic FWI.
Viscoacoustic and Viscoelastic Full Wavefield Inversion
FWI is a partial-differential-equation-constrained optimization method which iteratively minimizes a norm of the misfit between measured and computed wavefields. Seismic FWI involves multiple iterations, and a single iteration can involve the following computations: (1) solution of the forward equations, (2) solution of the adjoint equations, and (3) convolutions of these forward and adjoint solutions to yield a gradient of the cost function. Note that for second-order optimization methods, such as Gauss-Newton, the (4) solution of the perturbed forward equations is also required. A more robust mathematical justification for this case can be found, for example, in U.S. Patent Publication 20130238246, the entire content of which is hereby incorporated by reference.
For viscoacoustic and viscoelastic seismic FWI, the cost of each simulation is exacerbated due to the fact that the relaxation mechanisms needed to accurately model the quality factor's frequency-dependent profile introduce even more variables and equations. Furthermore, the number of simulations that must be computed is proportional to the total number of iterations in the inversion, which is typically on the order of hundreds to thousands. Nevertheless, the benefit of inferring the attenuation properties of the subsurface using this method is expected to outweigh the cost, and development of algorithms and workflows that lead to faster turnaround times is a key step towards making this technology feasible for field-scale data, allowing users to solve larger scale problems faster.
A common iterative inversion method used in geophysics is cost function optimization. Cost function optimization involves iterative minimization or maximization of the value of a cost function (θ) with respect to the model θ. The cost function, also referred to as the objective function, is a measure of the misfit between the simulated and observed data. The simulations (simulated data) are conducted by first discretizing the physics governing propagation of the source signal in a medium with an appropriate numerical method, such as the finite difference or finite element method, and computing the numerical solutions on a computer using the current geophysical properties model.
The following summarizes a local cost function optimization procedure for FWI: (1) select a starting model; (2) compute a search direction S(θ); and (3) search for an updated model that is a perturbation of the model in the search direction.
The cost function optimization procedure is iterated by using the new updated model as the starting model for finding another search direction, which will then be used to perturb the model in order to better explain the observed data. The process continues until an updated model is found that satisfactorily explains the observed data. Commonly used local cost function optimization methods include gradient search, conjugate gradients, quasi-Newton, Gauss-Newton and Newton's method.
Local cost function optimization of seismic data in the acoustic approximation is a common geophysical inversion task, and is generally illustrative of other types of geophysical inversion. When inverting seismic data in the acoustic approximation, the cost function can be written as:
                                          ℱ            ⁡                          (              θ              )                                =                                    1              2                        ⁢                                          ∑                                  g                  =                  1                                                  N                  g                                            ⁢                                                          ⁢                                                ∑                                      r                    =                    1                                                        N                    r                                                  ⁢                                                                  ⁢                                                      ∑                                          t                      =                      1                                                              N                      t                                                        ⁢                                                                          ⁢                                      W                    ⁡                                          (                                                                                                    ψ                            calc                                                    ⁡                                                      (                                                          θ                              ,                              r                              ,                              t                              ,                                                              w                                g                                                                                      )                                                                          -                                                                              ψ                            obs                                                    ⁡                                                      (                                                          r                              ,                              t                              ,                                                              w                                g                                                                                      )                                                                                              )                                                                                                          ,                            (        6        )                            where        (θ)=cost function,        θ=vector of N parameters, (θ1, θ2, . . . θN) describing the subsurface model,        g=gather index,        wg=source function for gather g which is a function of spatial coordinates and time, for        a point source this is a delta function of the spatial coordinates,        Ng=number of gathers,        r=receiver index within gather,        Nr=number of receivers in a gather,        t=time sample index within a trace,        Nt=number of time samples,        W=norm function (minimization function, e.g. for least squares function (x)=x2),        ψcalc calculated seismic data from the model θ,        ψobs measured seismic data (pressure, stress, velocities and/or acceleration).        
The gathers, data from a number of sensors that share a common geometry, can be any type of gather (common midpoint, common source, common offset, common receiver, etc.) that can be simulated in one run of a seismic forward modeling program. Usually the gathers correspond to a seismic shot, although the shots can be more general than point sources. For point sources, the gather index g corresponds to the location of individual point sources. This generalized source data, ψobs, can either be acquired in the field or can be synthesized from data acquired using point sources. The calculated data ψcalc on the other hand can usually be computed directly by using a generalized source function when forward modeling.
FWI attempts to update the discretized model θ such that (θ) is a minimum. This can be accomplished by local cost function optimization which updates the given model θ(k) as follows:θ(i+1)=θ(i)+γ(i)S(θ(i)),  (7)where i is the iteration number, γ is the scalar step size of the model update, and S(θ) is the search direction. For steepest descent, S(θ)=−∇θ(θ), which is the negative of the gradient of the misfit function taken with respect to the model parameters. In this case, the model perturbations, or the values by which the model is updated, are calculated by multiplication of the gradient of the objective function with a step length γ, which must be repeatedly calculated. For second-order optimization techniques, the gradient is scaled by the Hessian (second-order derivatives of objective function with respect to the model parameters). The computation of ∇θ(θ) requires computation of the derivative of (θ) with respect to each of the N model parameters. N is usually very large in geophysical problems (more than one million), and this computation can be extremely time consuming if it has to be performed for each individual model parameter. Fortunately, the adjoint method can be used to efficiently perform this computation for all model parameters at once (Tarantola, 1984). While computation of the gradients using the adjoint method is efficient relative to other methods, it is still very costly for viscoacoustic and viscoelastic FWI.
Conventional viscoacoustic and viscoelastic full-waveform inversion methods in the time domain compute the gradient of the memory variables to integrate the adjoint equations in time. Earth models taking attenuation into account have been presented both in the frequency and time domains (Ursin and Toverud, 2002). The main focus of frequency-domain methods has been establishing a relationship between attenuation and medium velocities using complex number properties under causality principals, which make them straightforward to directly apply for FWI (Hak and Mulder, 2010).